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The advantages of performing Langevin Dynamics in extended systems are discussed. A simple 
Langevin Dynamics scheme for producing the canonical ensemble is reviewed, and is then extended 
to the Hoover ensemble. We show that the resulting equations of motion generate the isobaric- 
isothermal ensemble. The Parrinello-Rahman ensemble is then discussed and we show that despite 
the presence of intrinsic probability gradients in this system, a Langevin Dynamics approach samples 
the extended phase space in the correct fashion. The implementation of these methods in the ab- 
initio plane wave density functional theory (DFT) code CASTEP [M. D. Segall, P. L. D. Lindan, 
M. J. Probert, C. J. Pickard, P. J. Hasnip, S. J. Clarke and M.C. Payne, J. Phys.: Cond. Matt. 14 
(11), 2717 (2003)] is demonstrated. 
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I. INTRODUCTION 

The method of molecular dynamics is a well established 
tool with wide ranging applications. Of particular advan- 
tage is the ability to obtain ensemble averages of statisti- 
cal quantities (in ergodic systems) and detailed trajectory 
information within the same computational framework. 

Many of the early attempts to extend molecular dy- 
namics from the micro-canonical (NVE) to the canonical 
(NVT) ensemble used some form of stochastic dynamics. 
A random component to the particle dynamics stimulates 
ergodicity and can reduce correlations for increased sam- 
pling efficiency at the expense of accuracy in short-term 
dynamics. A simple stochastic dynamics scheme based 
on the Langevin equation will be reviewed in section ITTl 

A more widely adopted scheme for generating NVT 
ensemble trajectories is the Nose thermostat |l|, U • Here 
time is rescaled according to an extra 'extended system' 
variable. This is introduced to the Hamiltonian in a man- 
ner chosen to reproduce the correct NVT partition func- 
tion and hence simulate coupling to a heat bath. Use 
of this method requires either a non-uniform sampling in 
the scaled time variable or a non-canonical transform 
to un-scaled time 0- Recent work on non-Hamiltonian 
statistical mechanics has shown that this is not a major 
concern for the trajectory of the particle subsystem 
Other work has reformulated the Nose scheme to remove 
the need for a non-canonical transform ||| altogether. 

Of greater concern is the problem of ergodicity. In 
order to generate correct ensemble averages from MD 
trajectories, the simulated system must be ergodic. The 
simple deterministic nature of the Nose scheme requires 
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that the chaotic behaviour needed to stimulate ergodicity 
is provided by the particle sub-system itself. For certain 
systems, in particular those containing harmonic forces, 
this requirement is known not to be satisfied Q leading 
to incorrect phase-space trajectories and hence poor en- 
semble averages. It is now common practise to overcome 
this problem by the use of a chain of Nose-style ther- 
mostats 0. Although guaranteed to correctly sample 
the isothermal ensemble, this scheme is still determinis- 
tic and exhibit longer correlation times than stochastic 
methods in some circumstances. Note that recent work 
has generalised the Nose scheme to include more chaotic 
terms in the extended Lagrangian 0, Q which does much 
to eliminate these issues. 

Extended Hamiltonian methods are also available to 
simulate the coupling of the particle system to a pressure 
piston. This idea was introduced by Andersen |^ in a 
scheme where particle positions and momenta are scaled 
according to the cell volume. Potential and kinetic en- 
ergy terms for the volume are added to the Lagrangian 
resulting in an expansion or contraction of the system to 
regulate pressure with the required fluctuations. A mod- 
ification of the Andersen scheme in which the extended 
system variables are the strain rate Msilon and its con- 
jugate, was introduced by Hoover [13 ■ This tends to 
be more robust in practise. A scheme in which the cell 
motion is anisotropic with each cell vector evolving in- 
dependently was introduced by Parrinello and Rahman 

Sia. 

These constant pressure schemes have been coupled to 
Nose style thermostats resulting in schemes for sampling 
the isobaric-isothermal (NPT) ensemble 

Kolb and Dunweg have shown that an effective 
method for sampling the isobaric-isothermal ensemble 
can be obtained by performing Langevin dynamics in 
the extended Andersen system. This has the poten- 
tial advantage of shorter correlation times and guaran- 
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teed ergodicity in the NPT ensemble. In this paper, we 
show that Langevin dynamics can be conducted within 
a Hoover-style extended system to correctly sample the 
isothermal-isobaric ensemble. The extension of these 
ideas to Langevin dynamics in a Parrinello- Rahman sys- 
tem will also be discussed. The resulting equations of 
motion are integrated using an evolution algorithm ob- 
tained from Louiville time evolution operators 16] which 
is symplectic in the limit of zero friction. 

Criteria for selecting coupling parameters for the ther- 
mostatic/barostatic processes will also be discussed and 
applied to various example simulations. 



II. LANGEVIN DYNAMICS 

A simple but effective method of performing Langevin 
dynamics simulations for the isothermal ensemble uses 
the following equations of motion in the usual notation 

U = pi/m,, (la) 
Pi = fi - jpi + Ri. (lb) 

Here 7 is a friction co-efficient representing viscous 
damping due to fictitious 'heat bath' particles. The ran- 
dom force Ri represents the effect of collisions with these 
particles. Following Chandrasekhar we assume that 
the timescale of the coUisional heat bath process is very 
much smaller than the ionic motions of interest, and 
hence 



i=l 

N 

= 7 51 ^p. • Ip^P + mksTVp^p] . (4) 

This has the canonical phase space probability density 
function as a stationary solution, and hence the 

method can be used to sample the isothermal ensemble 
via a fluctuation- dissipation processfigj. 

The choice of the parameter 7 is a compromise between 
statistical sampling efficiency and preservation of accu- 
racy in short-term dynamics. However, in the case where 
the process approximated by the stochastic components 
can be simulated by another method, it is possible to de- 
termine an optimal value of 7 numerically. Guidelines for 
the choice of 7 will be discussed in detail later. 

An effective Verlet-style integrator can be obtained for 
this scheme by applying the time-evolution operator for- 
mulation of Tuckerman et al in the limit 7 — > 0. The 
lack of a Louiville operator for the stochastic components 
of the dynamics requires that these are then included into 
the particle forces via the following substitutions: 



f, (i) ^ f, (i) - 7p, (0 + R, (i) 
fi (t + M) ii{t + At) - 7Pj (< + A<) + Rj (i + At) 



(R, (t) R.,it'))^5{t-t'). 



(2) 



Furthermore, we assume that a great many collisions 
with heat bath particles occur over a MD time-step and 
we can therefore expect that R^ will be distributed in a 
Gaussian fashion in accordance with the central limit the- 
orem. The Stokes-Einstein relation for the diffusion co- 
efficient can then be used to show that the average value 
of Rj; over a time-step (in thermal equilibrium) should be 
a random deviate drawn from a Gaussian distribution of 
zero mean and unit variance scaled by 



At 



(3) 



Note that this choice of R^ limits this simple scheme, 
and the constant pressure schemes that follow, to equilib- 
rium simulations only. Certain deterministic thermostats 
can however be utilised in non-equilibrium simulations 

El. 

It can be shown that eauations llal and ll bl are equivalent 
to the following Fokker-Planck equation for the phase 
space probability density p (r^,p^) 



In the NVT case described here, this leads to standard 
Velocity- Verlet, modified with the above substitutions. 



III. LANGEVIN-HOOVER DYNAMICS 

In this section, we show that the simple NVT scheme 
above can be extended to perform Langevin dynamics in 
a Hoover-style extended system, and that this results in 
correct sampling of the isobaric-isothermal ensemble. 



A. Equations of Motion 

The equations of motion we propose are shown below. 
These incorporate the constant pressure component of 
the Hoover equations as corrected by Martyna et al rTsj . 
In the form shown below, the deterministic equations for 
the evolution of both the particle and barostat velocity 
have been converted to Langevin equations in d dimen- 
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sions with difTerent friction constants 



V 

P.psilo7i 



Pi , P.psilon 

^ TT/ ' 

nii W 

-Vr,.$- (l- 



Nf 



p.psilon 

w 



dVp. psilon 

/w 



(5a) 

Pi - 7Pi + R(5b) 

(5c) 



dV (X - Pe^t) + TT — - IpP-psilon +(S^) 



■' ? — 1 



These equations introduce the volume V as a dynam- 
ical (barostat) variable. The corresponding momentum 
variable p.psUon is the strain rate .psilon multiplied by a 
fictitious mass W. Rp is a stochastic 'force' which acts 
on the barostat. The use of a Langevin equation for the 
barostat as well as the particles may have possible equi- 
libration benefits, but we shall see in the analysis that 
follows that it is not required for canonical sampling at 
equilibrium. 

The scalar X is given by 



where H. is the Hamiltonian of the particle sub-system. 
Note that H' itself is not a Hamiltonian. A rigorous 
phase space analysis is therefore best performed using 
techniques in non-Hamiltonian statistical _mechanics as 
introduced by Tuckerman and co-workers 
B. Justification 



To show that equations [53 to I5dl correctly sample 
the isobaric-isothermal ensemble we first identify the un- 
thermostatted phase-space probability density for the ex- 
tended system of particles, 



Pnph' (r^,P^,P.ps«ion, V) oc 



5[H'{t)-H'm 

^NPH' (r^,P^,P.psiion, V]' 
(10) 

The phase space density for this system coupled to a 
heat bath at constant temperature should therefore be 



X 



1 

dV 



N N 

Pi Pi , n 



,i=l 



i=l 



dV 



V). (6) 



The distinction between X and the instantaneous pres- 
sure V is important. The value of V when calculated 
using the virial equation should include the white noise 
contributions from the 'Langevin thermostat' whereas X 
does not. Use of V in eauation l5dl leads to non-canonical 
temperature and volume fluctuations. 

The values of are drawn from the same distribution 
as the NVT case. Values of Rp are drawn from a Gaussian 
distribution of zero mean and unit variance scaled by 



2kBTW-fp 



At 



In the un-thermostatted limit (7 0,Jp 
tions I5al to 15 bl obev the Louiville theorem 



(7) 

0) equa- 



Pnpt iy iP I P.psiionj 



X exp 



VInPT (r^, P^ .P.psilon, V) 

n{v^,Y>'')+PV + p%,aonl'^W 



(11) 



This is the probability density function for the cor- 
rectly thermostatted extended particle plus barostat 
phase space. Integration over the barostat momentum 
P.psilon yields a constant, and hence 



Pnpt (i"^,P^, V) oc 



1 



X exp 



n^PT (r^,P^,V) 



(12) 



^ TV N a o 

op \ - . _ \ - . _ . op op 

-qI + 2^^i^ yiP + 2^Pi^ PiP + P.ps^lon — + V— = 



and conserve the quantity 



' 9pe dV 



(8) 



i/' = H(r^,p^)+PV+p2 (9) 



which is the correct probability density function for the 
isobaric-isothermal particle subsystem. 

Following Stratonovich [23| we now construct the fol- 
lowing Fokker-Planck equation for the phase space den- 
sity p resulting from equations I5al to l5dl The compress- 
ibility of equations |53 to I5dl is zero and hence the Jaco- 
bian of the resulting co-ordinate transform is the identity 
matrix and we need not include it: 
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dp 

dt 



N 

E 



, P.psilon 



rfV(x-Pe.t) + l^E- 



- 1 1 

dp 



dp 



Nf 

psilon ^ ^ 



p.psilon 



dp.ps 



P.psilonP + WksT 



d 



dp.ps 



-P 



N 
i=l 



(13) 



where p = p (r^ , , p.psiion , V) . 

In order for the Langevin-Hoover scheme to correctly 
sample the isobaric-isothermal ensemble, equation 1111 
must be a solution of ^1 As the un-thermostatted sys- 
tem is Louivillian, we expect the LHS of eauation ll3l to 
be identically zero which is easily confirmed. Our use of 
the Stokes-Einstein relation in balancing the stochastic 
part of the Langevin dynamics with the dissipative fric- 
tion term ensures that the RHS is also zero for any phase 
space probability density function representing equilib- 
rium with a heat bath. This is also easily confirmed in 
the case where p is given by equation II II and hence the 
proposed equations of motion are expected to sample the 
isobaric-isothermal ensemble in the correct fashion. 

A complete justification that the Langevin-Hoover 
scheme samples the required ensemble must show that 
the desired temperature and pressure are maintained 
with canonical fluctuations about the mean of each. This 
behaviour depends on the choice of the parameters W , 7 
and 7p, and will be covered in detail in section Fill Dl 



C. Numerical Integration 

In this section we follow the method discussed above 
for obtaining integration algorithms for systems obey- 
ing Langevin equations. The Louiville operator for the 
isotropic Langevin NPT equations of motion in the limit 
of zero friction is 



d . d -,9 
+ + -Psilon— 
dr dp de 

iL,- + iLp + iL. psilon + i^Pe 



iL = i^ + t>^ y -psilon— + p.psiionTT- 

dpe 



(14) 



The following Trotter factorisation of the resulting 
time-step evolution operators was found to be the most 
convenient: 



1. yt+h^t = V* + 

t+^At _ ( At • 
^- P.vsilon — Pe + 2 P-Psilon 



.psilon 



.psilon 

r* n* V*+^^* n* ., 

J J f J I " I F.pstlc 



+ f p. 



3. p 

4. r*+'^* = r* -I- Atr 



I Pi : P .psilon 



t t+iAt t+iAt 
^2 ' Pi ' P .psilon 



0- Pi — Pj + — Pi 



't+At t + ^At t+^At 
i ' r'i ' P.psilon 



6. p 



t+At 

'silon 



t+iAt 
-Pe 



~t~ 2 P.psilon 



t+At „t+At yt+^At 

i ' r'i ' 'i^. psilon 



Making the appropriate substitution for the particle 
forces at each time-step we now include the Langevin 
buffeting and damping terms and denote these as extra 
dependences of the time derivatives. 



V\p* 



.psilon 



' ^ .psilon 
2 P.psilon 



o t+|At f , At • 

3- P, ' = P' + —Pi 



^iiHii" ^ P .psilon^ IPP .psilonT ''p 



n,Pi,P 



psilc 



,7P*,R* 



t+At „ j.t 



4. 

5. P^^* 



Atr-i 

t+iAt 



J f+iAt t-l-iAt 

'^iiPi 'P.psilon 



2 P» 



t+At t+^At t+iAt , 
'■I ' f I ' f. psilon I i^i ' '■i 



^t+At -rit+At 



giLAt ^ giL.p„,„„At/2gi-Lp p^;,^^At/2giLpAt/2giL,At^ 

giLpAt/2gi-Lpp^.,„„At/2giLp^,,„„At/2 (--^g^ 

which leads to the following integration algorithm. 



6 D*+'^* 
' p.psilon 



t+^At 



Pe 



2 P.psilon 



t+At „t+At yt+iAt „t+At „t+At nt+At 

i ' i^i ' 'P.psilon' ipP. psilon' ^p 
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FIG. 1: Value of the conserved quantity in the Hoover system 
for a Langevin dynamics simulation, and a similar simulation 
in the zero friction limit. Interactions were modelled with the 
Lennard- Jones potential for argon at, 67.6 MPa. In the case 
of the Langevin dynamics simulation, a set temperature of 
20 K was used. An equilibration time of 10,000 At was used 
prior to sampling. The drift in the zero friction value is less 
than 2 x 10~^ Hartree over the 100,000 time-step run with 
At = 2.4 fs. 



yt+iAi t+At 
' ^ .ps%lon 



To demonstrate the stability of this algorithm we con- 
sider the limit of vanishing friction in which 7^0. In 
this limit the system is Louivillian and the quantity 



N 



i=l 



+ pLzio„/2W^ + ^-tV (16) 



is conserved. Figure 1 shows the value of H' for a 
Lennard- Jones system both in the zero friction limit, and 
for a Langevin- Hoover simulation. With the introduction 
of finite friction we expect deviations from H' with no 
long term drift as indicated. 



D. Choice of Parameters 

The isotropic scheme uses three parameters which 
must be chosen carefully. These are the barostat 'mass' 
W , the particle friction coefficient 7 and the barostat 
friction coefficient 7^. 

It has been previously noted 21] that for a Hoover 
barostat, the fictitious mass should be chosen according 
to 



W = dNhT/ul 



(17) 



where lo^ is the frequency of the required volume fluctu- 
ations. The choice of this uJb depends on the timescale 
on which the particle motions of interest occur. A use- 
ful aide when choosing this frequency is the NVT tem- 
perature spectrum (TS). The character of this spectrum 
should not be significantly disrupted by the addition of 
a barostat. In practise this means choosing Wf, to be less 
than the smallest frequency component in the TS (figure 
2) . Care must also be taken that LOb is sufficiently close to 
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FIG. 2: Choice of uJb for solid Lennard- Jones argon at 50 K 
20 MPa. A thermostat frequency u]p = 2-K'y — 6.6 x 10~* rad 
fs~^ is used in both cases. In the left hand case, ujf, is chosen 
to be less than the lowest frequency which appears in the 
NVT temperature spectrum. The corresponding temperature 
and volume sample distributions match the exact canonical 
(solid line) case exactly. In the right hand case the choice of 
UJb has disrupted the natural temperature spectrum resulting 
in non-canonical fluctuations. 



these components that de-coupling of the barostat from 
the particle subsystem is avoided. 

To determine an appropriate value for the particle fric- 
tion coefficient 7 we calculate the memory function for 
the system we wish to study. The memory function ^ for 
the velocity autocorrelation function tp is defined by 



^{t~T)'4> (r) dr. 



(18) 



The generalised Langevin equation is 



p, = F,- f Ut-r) p^ (r) dr + R„ (19) 

JQ 

which in the extended particle + barostat phase space 
becomes 



psilon 



Nf J W 



■ / Ht-r)p^{T)+■R,. (20) 
^0 



6 



Comparing to eauation lSbl we can see that the stochas- 
tic component of the dynamics is generated within the 
approximation 

^ = l5{t). (21) 

To retain consistency with this approximation, the 
value of gamma should be equal to 



7 = ^0= / Ut{r)dT, (22) 
Jo 

where ^act is the actual memory function of the system we 
wish to simulate. Calculation of the optimal 7 therefore 
requires prior knowledge of the memory function which 
can be obtained from an NVE simulation or from other 
NVT/NPT techniques. As the stochastic component of 
the dynamics has now been chosen to be characteristic of 
the particle motions themselves, the 'heat bath' is now 
representative of the effect of the bulk on our simulated 
sample. This is particularly appropriate for simulations 
using periodic boundary conditions, where the only phys- 
ical objects with which the simulated particles exchange 
heat, are other particles in neighbouring cells. 

A useful method of calculating ^{t) from trajectory in- 
formation via an autoregressive model has been proposed 
by Kneller and Hinsen 22]. We will use this method to 
calculate appropriate value of 7 for various simple sys- 
tems in section fV Al 

The cost of computing a memory function in advance 
of a Langevin dynamics simulation can sometimes be 



prohibitive and it is often convenient to adopt the con- 
servative rule of thumb that the thermostat frequency 
LUp = 27r/7 should be a few tens of times smaller than 
the smallest characteristic frequency of the system to be 
studied. This may result in less than optimum sampling 
efficiency, but will generate the NPT ensemble without 
disrupting the particle motions of interest. 

The choice of the 'piston thermostat' parameter jp is 
of less importance. The inclusion of the damping and 
buffeting terms in eauation lSdI is not necessary for repro- 
duction of the NPT partition function, however a non- 
zero 7p may equilibration times. Furthermore the effect 
on the particle velocities is small in equilibrium. We gen- 
erally require that In j^p be ten times smaller than the 
frequency associated with the barostat variable. It is in 
principle possible to calculate an autocorrelation func- 
tion for p.psiion and hence a memory function and opti- 
mal barostatic friction coefficient. This is not useful in 
practise. 

IV. LANGEVIN-PARRINELLO-RAHMAN 
DYNAMICS 

A. Equations of Motion 

An implementation of Langevin Dynamics in a fully 
flexible simulations cell follows from the Nose-Hoover 
thermostatted Parrinello-Rahman scheme of Martyna, 
Tobias and Klein after removing the Nose-Hoover 
chains and converting to Langevin equations in the mo- 
menta: 



h 



Pi 



-Vr..<f (r^,h)- 



Wa 



1 \ Tr[pg] 



NfJ Wg 



Pi - 7Pi + R-i 



Pg = V(X-P,,iI) 



1 " n2 

_L Pi 

iV/ 'H' rai 



I - IpPg + R-p 



(23a) 
(23b) 
(23c) 

(23d) 



where the tensor X is given by 



V 



N 



,1=1 

.N 



a, (3 



(24) 



As before, we draw the random forces from the same 
distribution as in the NVT case. Each component of the 
barostat buffeting tensor Rp is drawn from a Gaussian 
distribution on unit mean and zero variance scaled by 



2kBWg-fp 
At ■ 



(25) 



7 



where Wg is set by the barostatic frequency ujh according 
to 



Wg = {Nf + d)kBT/dLol 



(26) 



These equations can be evolved using an analagous in- 
tegration scheme to that described for the isotropic case 
in section Fill (]l Approprate frequencies for the thermo- 
stat and barostat are also chosen in the same way as for 
the isotropic scheme. 



B. Analysis 
Equations to |23d| 

conserve the quantity 



H' = n 



(r^,P^) 



2^Tr[p,pJ] +Pe.tdet[h] (27) 



in the limit of zero friction (the pure Parrinello- Rahman 
system). As in the Hoover case, this is not a true Hamil- 
tonian. The system of equations cannot be derived from 
equation 1271 via Hamiltons equations. The phase space 
analysis therefore requires the tools of Tuckerman et al 
and merits more attention than in the Hoover case. 



The compressibility k of equations I23al to I23dl in the 
zero friction limit determines the uniformity of the back- 
ground system in which we wish to perform Langevin 
Dynamics. For this system 



N 



N 



d d d (h) 



d d 



(28) 



which can be simplified to 



K = (d - l)Tr [pg] /W^ - dTr [pg] /VK + dTr [pg] /VF + 
= (d - l)Tr [pg] /M^. (29) 

We then identify the scalar strain rate .psilon as 
Tr [pg] /Wd, and can therefore construct the Jacobian of 
the co-ordinate transform which takes the system from 
t to t' = t as 



J(t;0) = exp( / d{d — \) .psilon dt' 



(30) 



and hence the phase space metric is 



V5M) 



d-l 
ref 



1 



V 



d-1 



oc 



det (h) 



d-l ■ 



(31) 



where Vref is the volume of the cell to which the strain 
is referenced. We therefore expect the correctly ther- 
mostatted probability distribution in the extended phase 
space to be 



p (r^, p^, h, Pg) oc exp [-H'/ksT] det [h] 



l-d 



(32) 



Integration of this over the d^ components of the strain 
momentum tensor Pg again yields a constant, yielding 



.N „N 



X exp [-{H + Pext det [hD/fcsT] det [h]'"'' , 

(33) 

which has previously been identified |l3| as the correct 
phase space distribution for the isobaric-isothermal en- 
semble with a fully flexible cell. We can therefore con- 
clude that the Parrinello-Rahman Langevin dynamics 
scheme will be capable of correctly sampling this ensem- 
ble if equation is a solution of the following Fokker- 
Planck equation obtained from equations I23al to I23dl 



dp 



N 



N 



TO,; Wn 



f 

N 



1 \ Tr[pg]_ 



NfJ Wg 



(Xq^ - PextSap) det [h] -f -i- — 



Nf nii 



dp 



dp 



7 ^ Vp^ • [p^p + mkBTVp^p] + 7p X! 

1=1 a,/3 



diPg)al3 \Wg J d{h)af3 

dp 



{Pg)a0P + WgksT 



(34) 



r 



Analysis of the LHS yields zero. This is equivalent to generalised Louiville theorem in the limit of zero friction, 
the statement that the extended phase space obeys the The background system in which the Langevin dynam- 
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Time (ps) Tijuc (ps) 



FIG. 3: NVE (solid line) and NVT Nose-Hoover (dashed line) 
VACF and memory function for Lennard- Jones argon at 80 K 
3.36 GPa. These are computed from 256 atom trajectories us- 
ing the autoregressive model of Kneller and Hinsen (order 
150) as implemented in the nMoldyn program |23| . 




Time (ps) Tijuc (ps) 

FIG. 4; VACF and memory function for 4000 (sohd-line) 
and 256 (dashed-line) Lennard- Jones argon NVT Nose- 
Hoover simulations at 80 K and equal densities equivalent to 
3.36 GPa. 



ics simulation is conducted therefore has a constant, but 
non-uniform probability density which conserves equa- 
tion E] 

This introduces the possibility that the Stokes-Einstein 
relation which we have used to balance the diffusion and 
friction terms in our Langevin equations may be invalid. 
When balancing these two terms (equation|3Jl we have im- 
plicitly assumed that the background probability density 
is flat, i.e. no background probability gradients. Fortu- 
nately all such gradients in the Parrinello-Rahman sys- 
tem are perpendicular to the particle momentum axes 
and hence the first term on the RHS of equation 1341 rep- 
resents balanced Langevin dynamics and is zero. 

The second term of the RHS is however not zero. The 
background probability gradients do effect the balance of 
diffusion and friction for the Langevin equation in the 
cell momentum. In the case where symmetrised pressure 
and cell buffetting tensors have been used to eliminate 
cell rotations, the imbalance is small and proportional 
to 7p. Purists may therefore which to remove the cell 
thermostat. In practise, the effect makes little difference 
in real simulations and so it is possible to use a finite 
7p to aid equilibration. Any benefits this introduces are 
however, likely to be small. 



1. Step 1 

First we obtain a memory function relevant to this 
state point by conducting an NVE simulation (energy 
0.105 Hartree, volume (19.05A)3). 

The resulting VACF and memory functions are shown 
in figure 3. We also plot the equivalent functions cal- 
culated from a Nose-Hoover NVT simulation to ensure 
our calculated value of 7 will be appropriate for a system 
coupled to a heat bath. As can be seen the two mem- 
ory functions are near identical in this case. Numerical 
integration under the memory function reveals an appro- 
priate value for 7 of 2.41 x 10~'^rad fs~^. 

Note that in order to ensure that the value of 7 cho- 
sen is truly representative of a bulk liquid system, we 
must ensure that the simulation used in its identification 
covers a length scale larger than that of any spacial cor- 
relations. For the purposes of this example we therefore 
also conduct a simulation using 4000 atoms and compute 
a memory function. The result is shown in figure 4 and 
indicates that our value of 7 is suitable. 



V. EXAMPLES 
A. Lennard- Jones Argon 

In this section we will follow through the s.ps of con- 
ducting an isotropic constant pressure Langevin dynam- 
ics simulation for high pressure liquid argon modelled 
using the familiar Lennard- Jones potential. In particular 
we will investigate the temperature and density fluctua- 
tions of this system at 80 K, 3.36 GPa. This corresponds 
to a total energy of approximately 0.105 Hartree for 256 
atoms occupying a cubic cell of dimension 19.05 A. 

We will also simulate a solid argon system using the 
Parrinello-Rahman based scheme and again investigate 
the quality of the fluctuations. 



2. Step 2 

We now conduct an NVT Langevin dynamics simula- 
tion using the optimally identified value of 7 as described 
in section^ This allows a temperature spectrum to be 
calculated, providing a criterion for choosing a barostatic 
timescale. A run of 1638400 time-s.ps was performed in 
this case to provide a suitably illustrative example, how- 
ever much shorter runs can be used in practise. As this 
system is liquid, we do not expect any significant features 
in the spectrum as confirmed in figure 5, and we therefore 
choose the barostat to operate in low frequency region be- 
low the influence of the thermostat i.e. = 2 x 10""* rad 
fs-i. 
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FIG. 5: NVT (solid line) and NPT (dashed line) Langevin dy- 
namics temperature Fourier transforms for the Lennard- Jones 
system at 80 K 3.36 GPa. The barostat frequency is chosen 
so as to only affect the low frequency components introduced 
by the thermostat, and not the higher frequency components 
characteristic of the liquid itself. 



FIG. 7: Evolution of simulation cell during the 256 atom cubic 
cell LPR simulation of argon at 20 K 67.6 MPa. 



of a very high quality. 
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FIG. 6: Temperature and volume sample distributions cal- 
culated from Langevin-Hoover simulation compared to exact 
canonical cases for the Lennard- Jones argon system at 80 K 
3.36 GPa. The compressibility required to calculate the ideal 
volume distribution was obtained from an earlier Nose-Hoover 
based NPT calculation. 



3. Step 3 

With suitable values of LOp and ojb identified, a reliable 
NPT simulation can be conducted. We have now also 
identified the timescales associated with the problem and 
can safely increase the calculation time-step from 1.2 to 
9.6 fs. Using these parameters the system is equilibrated 
for 50,000 At and then sampled every 10 time-s.ps for a 
further 500,000. The resulting temperature and volume 
distributions are shown in figure 6. 



4-. Solid Argon with Full Cell Fluctuations 

To demonstrate the ability of the Langevin-Parrinello- 
Rahman (LPR) scheme to produce equally accurate tem- 
perature and volume sample distributions, we have also 
conducted simulations at 20 K, 67.6 MPa where argon is 
a stable solid. Again memory functions were calculated 
from appropriate NVE and NVT simulations, and an ap- 
propriate frequency for the barostat was identified from 
the temperature spectrum as above. This process re- 
vealed suitable parameters of ujp = 1.33 x lO^^rad fs^^ 
and oj^ = 4.17 x 10~^ rad fs"^ 

These values were used in a simulation of a cubic cell of 
FCC argon containing 256 atoms. This system was equi- 
librated for 20,000 At and then sampled every 10 At for 
a further 750,000. A time-step of 9.6 fs was used. The re- 
sulting evolution of cell vectors in shown in figure 7, and 
the calculated distribution of temperature and volume 
samples is shown in figure 8. As with the liquid simula- 
tions using the isotropic algorithm, the distributions are 



B. Ab-initio Silicon 

In this section we will use the LPR scheme for the more 
realistic application of simulating silicon. We will first 
use the semi-empirical potential of TersofF |23| to obtain 
an approximate memory function for crystalline silicon 
at room temperature and pressure. We will then use 
this memory function to choose parameters for a smaller 
silicon system using the LPR implementation in the ab- 
initio plane wave DPT code CASTEP 



1. Step 1 

Again we perform both NVE and NVT runs with 216 
atoms at the temperature and pressure of interest (using 
the Tersoff potential) in order to obtain a memory func- 
tion characteristic of the system in question. These are 
shown in figure 9. A time-step of At = 2.0 fs was used for 
a total run of 20,000 At sampling every 4 time-s.ps after 
an initial equilibration period of 10,000 At. The system 
was initialised with the expected equilibrium diamond 
structure. 

Integration over the memory function yields a suitable 
thermostat frequency of 3.05 x lO^^rad fs~^. This cor- 
responds to a thermostat period of 20.81 ps. 

For interest, we also compute a memory function for 
a Nose-Hoover NVT simulation at the same state point 
using CASTEP. This was calculated for an eight atom 
cell using a cut-off energy of 160 eV and a time-step of 
8 fs. The ultra-soft pseudo-potential method was used. 
The Brillouin zone was sampled at 4 k-points using the 
Monkhorst-Pack method, and the exchange and correla- 
tion functional was approximated using the LDA. The 
result is shown in figure 10. Although we expect this 
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FIG. 8: Temperature and volume sample distributions of 
solid FCC argon at 20 K 67.6 MPa calculated using the LPR 
scheme. Exact canonical distributions are shown as solid lines. 
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FIG, 9: NVE (solid line) and NVT (dashed line) VACF and 
memory function for the 216 silicon atom system modelled 
with the TersofF potential at room temperature and atmo- 
spheric pressure. 
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FIG, 10: Memory function and VACF calculated from the 
8-atom ab-initio simulation of silicon at room temperature 

and pressure. These are expected to exhibit severe finite size 
effects, but are similar to those obtained from the Tersoff 
potential in figure 9, 



memory function to be subject to significant finite-size 

effects, the thermostat period obtained from this mem- 
ory function is 16.95 ps, similar to that obtained from the 
larger system. 



2. Step 2 

An NVT Langevin dynamics simulation at the required 

temperature and pressure (using the Tersoff potential) of 
163,840 At sampled every 10 At yields the temperature 
spectrum shown in figure 11, For this example we use 
a value of uib set to 3.5 x 10~^rad fs~^ which is well 
separated from the dominant frequencies and a few times 
smaller than that of the thermostat. The temperature 
profile for a NPT run at this value of uJb is also shown in 
figure 11, indicating that the choice is appropriate. 



FIG. 11: NVT (sohd line) and NPT (dashed line) Langevin 
dynamics temperature Fourier transforms for the 216 atom 
silicon system modelled using the Tersoff potential. The baro- 
stat frequency has been chosen to avoid disruption of the nat- 
ural components and hence the two profiles are similar. 
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FIG. 12: Evolution of cell angles during the 8 atom ab-intio 
LPR simulation. 
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FIG. 13: Distribution of temperature and volume samples 
during the 8 atom ab-intio LPR simulation. The solid lines are 

the calculated canonical bulk values for comparison, although 
it should be noted that both distributions should differ from 
the bulk cases which are derived for large N. 



3. Step 3 

With appropriate thermostat and barostat timescales 
obtained from the classic;al potential, we now perform the 
ab-initio LPR simulation. Again we use 4 k-points and a 
plane wave cut-off of 160 eV. The system was initialised 
in the ideal BC8 structure (cell parameter 5.75 A) with 
random velocities, and equilibrated for 8,500 A t. The 
system was then sampled at a further 6,200 s.ps. This 
represents a relatively simple simulation, requiring mod- 
est CPU time. 

The resulting cell evolution is shown in figure 12. We 
also calculate the distribution of temperature and volume 
samples as shown in figure 13. For such a small simula- 
tion, we do not expect exact fits to canonical distribution 
functions, however the plots indicate convergence toward 
the ideal results as demonstrated for argon above. 



VI. DISCUSSION 

The results presented in the previous section have 
demonstrated that Langevin-Hoover and Langevin- 
Parrinello-Rahman schemes are capable of producing 
canonical fluctuations from which ensemble specific ther- 
modynamic quantities can be calculated. Clearly the 
choice of parameters is important. 

In this paper we have calculated optimal damping co- 
efficients (and hence thermostat frequencies Wp) for the 
stochastic component of the dynamics by computing a 
memory function for our example simulations. While 
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this has provided optimal samphng efficiency, it has also 
complicated the choice of the barostat frequency ojb. We 
have stated earlier that the barostat frequency should be 
smaller than, but close to that of the simulated parti- 
cles. This allows the particle motions to modulate the 
volume oscillation and hence generate canonical fluctu- 
ations. Too low a barostat frequency, and its motion 
becomes decoupled from the rest of the simulation, re- 
sulting in an undesirable slow harmonic oscillation of the 
volume. We also require the thermostat frequency Up to 
be significantly higher than the barostat frequency lui,. 
This prevents disruption of the thermostatic process by 
the cell motion. 

As we have chosen our Langevin friction co-efficient 7 
to be characteristic of the particle motions, Up will be 
close to the particle frequency lo. We also require Wb to 
be close to lu but well separated from ujp. The choice 
of an optimal friction co-efficient therefore introduces an 
apparent conflict in the choice of barostat parameter. We 
have shown above in section IlII Dl that this issue can be 
overcome by use of the temperature spectrum in carefully 
choosing the barostat parameter. 

It is important to note that this is not always the case 
in practise. The optimal friction co-efficient identified 
from the memory function represents the boundary be- 
tween performing molecular dynamics with a 'Langevin 
thermostat' (i.e. low tOp with similar correlation times 
and sampling efficiencies as a Nose-style thermostat 
for chaotic systems) and performing a 'true Langevin 
dynamics' simulation where the stochastic components 
dominate and are guided by the particle interaction forces 
(high ujp). For the purposes of sampling an isothermal 
ensemble, a small conservative choice of cUp will also pro- 
duce canonical temperature fluctuations, but on a longer 
timescale. The value calculated from the memory func- 
tion therefore represents the best compromise between 
efficiency and preservation of accuracy in short-term dy- 
namics. 

As a conservative rule of thumb, ujp can be set a few 
tens of times smaller than the particle frequency lu, and 
LUb can be chosen a few times smaller again. This allows 
canonical simulations in either the isotropic or fully flex- 
ible cell NPT ensemble to be generated with ease, albeit 
with a reduced sampling efficiency due to weaker modula- 
tion of the volume oscillation. A compromise must there- 
fore be reached between the computational expense of 
computing optimal parameters, and statistical efficiency 
in the simulation itself. 

This stochastic approach bypasses any concerns about 
possible 'hidden' conservation laws 4] in Nose style 
schemes which can incorrectly restrict trajectories and 
lead to incorrect ensemble averages. This effect in a 
simple harmonic potential is well documented (although 
still not entirely understood) and is generally thought to 



be ignorable in practical simulations. Often a chain of 
thermostats is used to break any unphysical conserva- 
tion laws, however this process takes place on a longer 
timescale than in a stochastic thermostatting scheme - 
an important consideration in ab-initio dynamics. 

Solids at low temperature are essentially harmonic, as 
are the forces used to 'connect' different realisations in 
imaginary time during path-integral molecular dynam- 
ics simulations. We therefore suggest that the isobaric- 
isothermal sampling schemes presented here are a desir- 
able alternative to traditional deterministic schemes for 
many applications. 



VII. CONCLUSIONS 

We have shown that performing Langevin dynam- 
ics within constant pressure extended systems is a 
valid method for simulating the equilibrium isobaric- 
isothermal ensemble. An analysis of the phase space in 
these simulations has been presented, which shows that 
the correct distribution of probability is generated, pro- 
vided the non-Hamiltonian nature of the extended system 
is accounted for. 

We have derived a suitable integration scheme, anal- 
ogous to the Velocity- Verlet scheme, which allows re- 
versible integration of trajectories in the zero friction 
limit, and eliminates long term drift of the conserved 
quantity in the extended system when the stochastic 
component is introduced. 

Suitable distributions of temperature and volume sam- 
ples can be achieved provided a little care is taken in 
choosing parameters. Methods for identification of opti- 
mal values of these parameters have been presented, We 
have also discussed less stringent criteria for choosing pa- 
rameters in the case where computing a memory function 
or an accurate NVT temperature spectrum is computa- 
tionally prohibitive. 

We consider this scheme to be a desirable alterna- 
tive to extended system NPT schemes with Nose-style 
thermostats, in particular for small or approximately 
harmonic systems, where deterministic schemes can suf- 
fer from unexpected coupling to the particle sub-system 
leading to incorrect trajectories, and consequent thermo- 
dynamics averages. 
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